FILE COPY — DO NOT TAKE 


NIST GCR 06-906 


Probabilistic Models for Directionless 
Wind Speeds in Hurricanes 


Mircea Grigoriu 
Cornell University 


NUT 


National Institute of Standards and Technology 
Technology Administration, U.S. Department of Commerce 


NIST GCR 06-906 


Probabilistic Models for Directionless 
Wind Speeds in Hurricanes 


Prepared for 

U.S. Department of Commerce 

Building and Fire Research Laboratory 
National Institute of Standards and Technology 
Gaithersburg, MD 20899-8611 


By 

Mircea Grigoriu 
Cornell University 
Ithaca, NY 14853 


December 2006 


U.S. Department of Commerce 


Carlos M. Gutierrez, Secretary 


Technology Administration 
Robert Cresanti, Under Secretary of Commerce for Technology 


National Institute of Standards and Technology 
William Jeffrey, Director 


cals Be | lobo abil 
pore kt | i; | ngsae bai 


ay ; 


wer “oe roe ~ or naa: 
“ > iin ay sei 


: + a 
7 f 


Contents 


1 Dim CWO CUT CON OU Bere a resol acocsbedadnvaceatecetiiaseess l 
2 Probabilistic models for hurricanes ....0....0..000..00..0ccccceee l 
2.1 Hurricane time model and design wind speeds ............................ I 
Zed Hurricane intensity. Gumbel distribution ..........00000 2 
2.3 Hurricane intensity. Shifted Gamma distribution ...................... 3 
2.4 Hurricane intensity. Reverse Weibull distribution. .................... 4 
25 Hurricane intensity. Generalized Pareto distribution ............... 5 
DDS ie Te EITC COC IOI INIOIMLGNNLS eeetece. eater tees erttes sence cesascccees ersessecccete catelocces. 6 

2.5.2. The method of probability weighted moments ............................e 6 

ce ME TDG TTVCCINOGLTO COED EEE AN AN os ate gnoc ns fa darn cis sencaddvscsesaace’s 6 

3 Monte .@arlo'alcorithmes 629 i ee eee 7 
INLAND PL by N SS CTT GUCCI) 1 ae ate ee en Pe 8 
5 CE OTICIUESIOTISG eet ee PTI ooo cacacesnenncnstiicssneet 9 
IRE CWT UREN Oe Are, es, Sag 17 i ee ne 9 
Appendix: MATLAB function hurr_nd_me.m .........000.000.0..... 10 


7 ? 
7 
_ 
J 
7 
re if ! Ti? « ef 


5 ‘ vaty 1 is a otett ar.cary ot 


-) .-toinl Caeahring 


iy . viak 


il 
é rf i 
" ~* sae THK } | } oy 
a ik 
TAS PARES 
7 


of 


iauleRe ; 


» 


Probabilistic models for directionless wind 
speeds in hurricanes 


Mircea Grigoriu 
Cornell University, Ithaca, NY 14853 


1 Introduction 


Extensive records of simulated hurricane wind speeds for 16 azimuth directions are 
available at, for example, the NIST web site http://www.nist.gov/wind. These records 
cover time periods of the order 2,000 years so that may not be sufficient for estimating wind 
speeds with average return periods of the order 1,000 years. Probabilistic models are needed 
to predict extreme wind speeds corresponding to such return periods. 

Our objectives are to (1) develop probabilistic models for hurricane wind speeds recorded 
irrespective of direction, referred to here as directionless wind speeds, (2) calibrate these mod- 
els to wind records, and (3) provide a Monte Carlo algorithm for generating data sets over 
long time periods that are consistent with a site statistics. Four models are considered for 
extreme wind speeds, the Gumbel, the shifted Gamma, the reverse Weibull, and the Pareto 
distributions. The study uses extreme wind speeds irrespective of direction, referred to as 
directionless wind speeds. 


2 Probabilistic models for hurricanes 


Hurricanes arrive at a site at random times and are characterized by random extreme 
wind speeds. Records can be used to estimate the mean number vy of hurricanes per year 
at a site. Records can also be used to construct the distribution F’ of extreme wind speeds 
under the assumption that wind speeds recorded in different hurricanes belong to the same 
statistical population, that is, they are independent samples of F’. 


2.1 Hurricane time model and design wind speeds 


Let T; < Ty <--- be a random sequence denoting the arrival times of hurricanes at a 
site, and let X,, Xo,... be hurricane wind speeds irrespective of direction. 

It is assumed that (1) the random sequences 7,7>,... and X,, X2,... are mutually 
independent, (2) the random variables X,, X2,... are independent and have the same dis- 
tribution F’, (3) the hurricane season begins on June 1 and ends November 30, and (3) the 
random variables 7;,7>—7},... are independent and follow an exponential distribution with 


mean 1/v, v > 0, so that T;,7>,... are points of a Poisson process N with intensity v. The 
Poisson process N is on only during the hurricane season. 

Consider a time interval [0,7], 7 > 0, and note that N(7) gives the random number of 
hurricanes in [0,7]. The distribution F’, of the largest extreme wind speed in [0,7], that is, 
the random variable max;<ij<n(7) {Xi}, is 


Por) exp [ —yr(1—- F(z))| CE) 


since the probability of n hurricanes occurring in [0,7] is P(N(r) = n) = (vT)” exp(—vr)/n! 
and the probability that the conditional random variable max;<;<y(r){Xi} | (N (T= n) does 
not exceed x is equal to F(x)”. 

Suppose we retain from the sequence X,, X2,... of wind speeds at a site only those 
values exceeding a threshold x, so that the resulting process describes hurricanes with wind 
speeds larger than x. The mean arrival rate of these hurricanes is v(z) = v (1 — F(z)), so 
that 1/v(x) gives the average time between consecutive wind speeds larger than x. Hence, 
the design wind speed that is exceeded on average once in r years, that is, the r-year wind 
speed, results by setting 1/v(x) equal to r and is 


1 
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x ( ~) (2) 


We consider four models for F’, the Gumbel, the shifted Gamma, the reverse Weibull, 
and the generalized Pareto distributions. The method of moments and other methods are 
used in the following two subsections to estimate the parameters of these models from wind 
records. The uncertainty in the estimates of the parameters of Ff and the corresponding 
r-year wind speeds is quantified in a subsequent section. 


2.2 Hurricane intensity. Gumbel distribution 


A random variable X is said to a Gumbel or extreme type 1 variable with parameters 
(a, u) if it has the distribution 


F(e) = exp | -exp(-a(e-w))], —00 < © <0, (3) 
and density 
f(2) =a exp| ~ a(x -w) ~exp(—a(e-w)], Soo aise (4) 


We use the notation X ~ EX1(a, u) to indicate that X is a Gumbel variable with parameters 
(a,u). The parameters (a, u) are denoted by (al_gumbel, u_gumbel) in the MATLAB code 
hurr_nd_mc.m. 
The mean yp and standard deviation o of X are related to the parameters (a, u) by 
0.577216 
== LL aS 
a 


=— (5) 
o V6 . 

These relationships with jz and o replaced by their estimates fi and 6 are used in hurr_nd_mc.m 

to obtain estimates of (a, u). 
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2.3 Hurricane intensity. Shifted Gamma distribution 


Let X be a Gamma random variable with parameters (k, A), k > 0, A > 0, and density 


\M gk-l eae 


LG) seer rene ©, (6) 


where ['(-) denotes the Gamma function. The distribution of X is 


rae / nine {| t% a rata SK) (7) 


and is given by, for example, the MATLAB function cdf(’Gamma’,z, k,1/A). The solution 
tp of F(x,) = p € [0,1], that is, the p-fractile of F can be obtained from the MATLAB 
function icdf(?Gamma’,p, k, 1/X). 

Consider the random variable X = a + X, where a is a real constant. Since FOX < 
x) = BOX <2“ - a), the density and distribution of X are given by Eqs. 6 and 7 with 7 —a 
in place of x, and are valid for z > a. The mean p, variance o”, skewness 73, and kurtosis 
va of X are 


w=at+k/r 
6 AO 
Y3 = 2/ Vi 
Ya = 3(14+ 2/k). (8) 
Suppose that the sequence of extreme wind speeds X 1, Xo,... follows a shifted Gamma 
distribution and that a record (2,...,%p) of this sequence is available. Our objective is 


to estimate the parameters (a,k, A) of this distribution. The maximum likelihood method 
and the method of moments are commonly used to calculate estimates (4, k, A) of (a, k, ). 
Maximum likelihood estimates have smaller variance that those obtained by the method of 
moments but can be unstable for some values of k ([3], Section 7.1). To avoid such difficulties, 
the method of moments is used to estimate (a,k,). Suppose that estimates (/1, 6, 73, V1) 
of the moments (1, 0, 3,774) of X1, X2,... have been calculated from the available record 
(%1,---,2n). Then k can be estimated from the expression of the coefficient of skewness or 
kurtosis and the estimates of these coefficients. For example, k= 4/ 73” if the relationship 
between k and 73 is used. Alternatively, k can be determined from the relationship between 
k and y4. The MATLAB code hurr_nd_mc.m uses the relationships 


i 4/ Ys" 

K = Vi/é (9) 
to find estimates of & and A, and calculates estimates a of the shift a by two options. The first 
option sets @ = minj<;j<p{z;}. The second option calculates @ from @ = ju—k/X. If the result- 
ing estimate of a is such that @ > minj<j<p{z;}, we set @ = minj<j<,{z;}. The parameters 


(a, k, X) are denoted in hurr_nd_mc.m by (shift1, lamq1, kq1) and (shift2, lamq2, kq2) 
for options 1 and 2, respectively. 


2.4 Hurricane intensity. Reverse Weibull distribution 


Let Y be a Weibull random variable with parameters a > 0, € € R, and c > QO, 


distribution P 
Te aGX Delees uf) r = 
F(y) = P| & YO (10) 
0 


Ee, 


fy) = = (way exp  - (HX) ]. Tees (11) 


Values of distribution Fy) and solutions of F(y,) = p € [0,1] can be calculated by, for 
example, the MATLAB functions F'(y) = cdf(’wbl’, y—€, a, c) and y,—€ = icdf(’whl’, p, a, c). 

Moments of any order of Y can be obtained from moments E [YY] =I(1+q/c) of the 
scaled random variable Y defined by Y =&+a Y ((3], Chapter 20). For example, the mean 
fly, Variance o;, and skewness 7,3 of Y are 


and density 


by = € +aT(1 +1/e) 
g, =a (T+ 2/c) -T(1 + 1/¢)’) 
oy) Uee s/o) > Sey Cie Zee iia ole 
V3 > [> he eee a oS Oe (12) 
(T(1 + 2/c) —-T(1 + 1/c)?) 


The properties of the random variable X £ —Y, called reverse Weibull variable, result 
from those of Y. For example, the distribution of X, that is, the probability P(X < x) = 
P(Y > —z), is equal to exp | — (%*)°*] for x < n = —€ and 1 for x > 7. 


Suppose that the sequence of extreme wind speeds Xj, Xo,... follows a reverse Weibull 
distribution and that a record (1,...,2n) of this sequence is available. Our objective is 
to estimate the parameters (a,7 = —€,c) of the reverse Weibull distribution. The method 


of moments, the method of maximum likelihood, the method of probability-weighted mo- 
ments, and other methods can be used to estimate the parameters of this distribution ((4], 
Chapter 22). Extensive numerical studies suggest that the method of moments delivers 
satisfactory estimators for the unknown parameters of F’, in contrast to, for example, the 
maximum likelihood method that can produce unstable estimators [6]. These features of the 
method of moments and its simplicity are the reasons for selecting the method for our anal- 
ysis. The following 3 step algorithm can be used to estimate the parameters (a,7 = —&,c) 
by the method of moments. 


Step 1. Construct the record (y; = —21,.--,Y%n = —Zn). Since the record (%j,..., Zn) 
is assumed to consist of independent samples of a reverse Weibull distribution with 
parameters (a,7 = —&,c), the record (y1,..-,Y%) consists of independent samples of a 


Weibull distribution with parameters (a, €,c). 


Step 2. Calculate estimates ji, oF , and 4,3 for the mean ju, variance oF , and skewness 
coefficient Yy,3 om the sample (41,98, ¥,)- For example, j,, =), 4.U;/te ne 
Dias (yi - [ty) /n are estimates of the central moments of order q > 2, oF = M2, and 


ae tg [G4 tore 


Step 3. Estimates the parameters (a,£,c) from Eq. 12. First, find an estimate ¢ for 
c from the last equality in Eq. 12 with 4,3 in place of 7,3. This nonlinear equation 
needs to be solved by iterations. Second, find an estimate @ for a from the second 
equality in Eq. 12 with (o;, c) replaced by (s; ¢). Third, find an estimate € for € from 
the first equality in Eq. 12 with (u,,a,c) replaced by (fuy, d, ¢). lie fags mini<i<n{y}, 
then set € = minj<j<n{y;} and calculate (4, @) from the first equalities in Eq. 12 with 
(fly, 67) in place of (uy, 07). 


The estimates (@, €,é) of the parameters (a,€,c) delivered by the above algorithm are 
denoted in hurr_nd_mc.m by (alw, xiw, cw). 


2.5 Hurricane intensity. Generalized Pareto distribution 


Let X be a (maximal) generalized Pareto variable with distribution 


re 1— (1—ka/a), li kajoe 0, kel Pa.0 (13) 
1 —exp(-—2/a), oe WE eth Ge kU 
and density 
AG) = (1/a) (1—k2/a)”, Wohi Ueks 0a > 0 (14) 
(ly G)exp (=s7/0) 0.60, a > 0, 


depending on the scale and shape parameters a and k, respectively. The range of the 
argument x of F is [0,co) and [0,a/k] for k < 0 and k > 0, respectively. We denote the 
distribution F’ in Eq. 13 by GPD(a,k). The parameters (a,k) of F and the mean yp and 
variance o” of this distribution are related by 


ay : (u?/o7 +1) 
1 


hes (u?/o? — 1). (15) 


The p-fractile rz, of F ~ GPD(a,k), that is, the solution of F'(x,) = p, has the expression 


a(1—(1—p)*)/k, k4#0 
p= ik ea ba: a) 


A notable property of generalized Pareto random variables is that, if X follow a 
GPD(a, k) distribution, then the conditional random variable (X —a) | (X > a) is GPD(a— 
ka,k) for any a € R. This invariance property simplify significantly the construction and 
analysis of peaks over threshold sequences associated with generalized Pareto series. 

Suppose that a record (21,...,@p) of independent values of X with distribution F' in 
Eq. 13 is available. Our objective is to estimate the parameters (a, k) of F from the record 
(%1,...,2%n). The maximum likelihood, probability-weighted moments, moments, and other 
methods can be used to estimate the parameters of F' ([1], Section 10.3 and 10.8). Extensive 
simulation studies indicate that estimates of (a, k) delivered by the method of moments is 
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generally reliable unless k < —0.2 and that the method of probability-weighted moments is 
adequate for k < 0 [2]. Other studies found the method of moments to be satisfactory for 
a broader range of values of k [5]. The MATLAB code hurr_nd_mc.m uses three methods 
for estimating the parameters (a,k) from data, the method of moments, the method of 
probability weighted moments, and the method of DEHAAN. 


2.5.1 The method of moments 


Estimates of the parameters (a, k) can be calculated from 


(17) 
where Z = )77_,2;/n and s? = )\, (a — Z)/n denote the sample mean and variance of 
(11,.-.-,%p), respectively. 

Suppose we construct from the original record (x,...,2n) a new record (y1,.--, Ym), 
m <n, consisting only of those readings x; exceeding a specified threshold a. The relation- 
ships in Eq. 17 with (m, {z; = yj — a}) in place of (n, {x;}) can be used to estimate the 
parameters of the Pareto variable Y — a. 


2.5.2 The method of probability weighted moments 


Let {2i-n} be the data set {z;} sorted in increasing order, that is, ty1.. < Lan < ++: < 
Inn. Set Din = (¢ —0.35/n), 2 ="1,..8, my and & —(1/n) 99 (et) 
estimates of a and k are given by ([1], Section 10.8.3) 


A tf 
oe 
Z—2t 
“ oa, 
INS 18 
%—2t aie 
Similar estimates can be constructed for a record (y,..-,Ym), m <n, consisting of 
values of (21,...,2n) above a threshold a. 
2.5.3 The method of DEHAAN 
Let (y1,---,;Ym), mM <n, a record consisting of values of (21,..., 2p) above a threshold 


a. Denote by {Yy;-m} the data set {y;} sorted in increasing order, that is, Yt:m < Yom <-+- < 
Ym:m, and consider the statistics 


1 Le 
pl) = —— [ In (Yj41:m) —Iln cheat to eae (19) 


for some q < n. The DEHAAN estimates of a and k are defined in the NIST web site 
http://www.nist.gov/wind and are given by 


; 1/2 
pA jane rnd sil anit 
tS (u(1)) |p?) 
a= ap! /p, (20) 


where p = 1 for k <0 and p=1/(1—&) fork >0. 

The estimates of the parameters (a, k) delivered by the methods of moments, probability 
weighted moments, and DEHAAN are denoted in hurr_nd_mc.m by (al_paretol, k_paretol1), 
(al_pareto2, k_pareto2), and (al_pareto3, k_pareto3), respectively. The threshold be- 
low which wind speed readings are disregarded is denoted by a_pareto for all estimation 
methods. 


3 Monte Carlo algorithm 
The hurricane hazard at a site during a time interval [0,7] is completely specified by 
— The arrival times 0 = To < Ti < T) <--- < Ty) of hurricanes and 
— The extreme speeds X1,X2,...,Xn(7) recorded during each hurricane. 


We have assumed that (1) the sequences {7;} and {X;} are independent of each other, 
(2) the arrival times {7;} define a Poisson process N with intensity v hurricanes/year, and 
(3) the random variables X,, X2,... are independent and follow either a reverse Weibull 
distribution or a shifted Gamma distribution. 

A two-step Monte Carlo algorithm has been developed to generate samples of extreme 
wind speeds recorded during hurricane at a site. The input to the algorithm consists of a 
reference time interval [0,7], the hurricane mean arrival rate v, the distribution shape F' of 
wind speeds X,, X2,..., and the parameters of this distribution. 


Step 1. Generate a sample 0 = To < Ti < Ty < --- < Tyiz of hurricane arrival 
times in [0,7]. As previously stated, it is assumed that hurricane can only occur 
from June 1 to November 30. Denote this time interval by 7,. The generation of 
hurricane arrival times can be based on the observation that the inter-arrival times 
Ty — Tr-1, k = 1,..., N(7) are independent exponential random variables with mean 
1/v. Hence, the time intervals 7; — T,_; are equal in distribution with the random 
variables — In(U;,)/v, where U; are independent random variables uniformly distributed 
in [0,1]. The MATLAB function rand can be used to generate samples of U;,. First, 
we generate a sample of 7; = —In(U;)/v. If this sample is larger than 7,, there will 
be no hurricane at the site in this sample during a yearly hurricane season so that 
N(ty) = 0. Otherwise, the sample of 7) gives the time of the first hurricane and 
N(ty) > 1. Second, we generate a sample of T2 — T; = —In(U2)/v. If the sample 
of T2 = (To — T,) + T; is larger than 7, then N(7,) = 1. Otherwise, N(7,) > 2 and 
we generate a sample of the following inter-arrival time. This process ends when an 


A 


arrival time exceeds 7, for the first time. A hurricane sample in [0,7] consists of the 
sequence of hurricanes generated in each hurricane season of [0,7]. Denote by N(r) 
the number of hurricanes in [0,7]. A similar procedure applies for the arrival times 
of hurricanes with wind speeds larger than a threshold ajp,;. The only difference is 
that the mean rate v is replaced with the mean rate v, equal to y scaled by the ratio 
#{hurricanes with speed > atnr}/#{all hurricanes}. 


Step 2. Generate a sample x = (21, 2%2,...,£n(r)) of the wind speeds X1, X2,..., X n(x) 
corresponding to the sample of 0 = Tp < Ti < Ty < --- < Tyi7) generated in the 
previous step. Let uw = (uj, U2,...,Un(r)) be independent samples of a uniformly 
distributed random variable in [0,1], which can be produced by, for example, the 
MATLAB function w = rand(1, N(7)). MATLAB functions have also been used to 
generate samples x of directionless hurricane wind speeds from uw under various as- 
sumptions on the distribution of wind speed. For example, the sample a is given by 
a at (€ + icdf(’wbl’, w, a, c)) for the reverse Weibull distribution. 


MATLAB function 


A MATLAB function hurr_nd_mc.m has been developed for estimating the parameters 


of the Gumbel, shifted Gamma, reverse Weibull, and generalized Pareto distribution from 
hurricane wind speed records. 


The input consists of: 
(1) A record at a specified milepost (see lines 50-56), 


(2) A range [cmin, cmax] for possible values of the tail parameter c in Eqs. 10-11 of 
the reverse Weibull distribution and the number nc of intervals to be considered in this 
range, 


(3) A threshold a_pareto for used for the generalized Pareto distribution, 


(4) A number nyr of years selected for Monte Carlo simulation and a seed nseed for 
generating hurricane wind speeds, and 


(5) A specified tail parameter cws for the reverse Weibull distribution. 
The output consists of: 
(1) A vector thurr with entries counting the number of hurricane in nyr years, 


(2) Vectors of simulated hurricane wind speeds: x_gumbel_mc corresponding to the 
Gumbel distribution; x_shgl_mc and x_shg2_mc corresponding to the shifted Gamma 
distribution under option 1 and option 2; x_rw_mc corresponding to the reverse Weibull 
distribution; and x_parl_mc, x_par2_mc, and x_par2_mc corresponding to the gener- 
alized Pareto distribution with parameters estimated by the methods of moments, 
probability weighted moments, and DEHAAN. | 


(3) Plots showing (2) histograms of the input wind record and Gumbel, shifted Gamma, 
reverse Weibull, and generalized Pareto densities fitted to this record, (77) wind speeds 
of average return periods in the range [50,1000] years predicted by the Gumbel, shifted 
Gamma, reverse Weibull, and generalized Pareto distributions fitted to the record using 
estimated parameters of these distributions and similar wind speeds derived directly 
from data, and (227) wind speeds of average return periods in the range [50,1000] years 
predicted by the reverse Weibull distribution with estimated and imposed tail param- 
eter and similar wind speeds derived directly from data. 


5 Conclusions 


A MATLAB code has been developed for estimating parameters of the hurricane wind 
speeds described by Gumbel, shifted Gamma, reverse Weibull, and generalized Pareto dis- 
tributions. Several estimation methods have been applied to calibrate these dsitributions to 
wind speeds records obtained from the NIST site hppt://www.nist.gov/wind. The result- 
ing models have been used to generate synthetic hurricane wind speeds and estimate wind 
speeds with return periods in the range [50,1000] years. The code produces numerous figures 
showing data and model statistics. : 


References 


[1] E. Castillo, A. S. Hadi, N. Balakrishnan, and J. M. Sarabia. Extreme Value and Related 
Models with Applications in Engineering and Science. Wiley Series in Probability and 
Statistics, Hoboken, New Jersey. 


[2} J. R. M. Hosking and J. R. Wallis. Parameter and quantile estimation for the generalized 
pareto dsitribution. Technometrics, 29(3):339-349, 1987. 


[3] N. Johnson and S. Kotz. Distributions in Statistics: Continuous Univariate Distributions- 
1. Houghton Mifflin Company, Boston, MA., 1970. 


[4] N. L. Johnson, S. Kotz, and N. Balakrishnan. Continuous Univariate Distributions, 
volume 2. John Wiley & Sons, Inc., New York, 1994. Second Edition. 


[5] E. S. Martins and J. R. Stedinger. Generalized maximum likelihood pareto-poisson 
estimators for partial duration series. Water Resources Research, 37(10):2551—2557, 2001. 


[6] Martins E. S. and J. R. Stedinger. Generalized maximum-likelihood generalized exteme- 
value quantile estimators for hydrological data. Water Resources Research, 36(3):737- 
744, 2000. 


Appendix. MATLAB function hurr_nd_mec.m 


function [thurr,x gumbel mc,x shgl me;x ‘shgZ2Z omc, x irw mC, 
x) parlime, x.pan2 Mc; X sparj Melee 
hurr _nd_mc(cmin, cmax,nc, Cws,a_pareto,nyr,nseed) 
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z The “Garectioniess wind dane ae tie (loupaney il 
% NEED FO MODEFY FOLLOWING INSTRUCTION 
TO) SENDING AMER EE: DESTRED MILEPOST # 


load miieposte6d 
% Poa) mil 
q=matrix; 
nr=length(q(:,1)); 

nu=mean_ rate; * mu = the average number of ee 


, ; : ; ? 
aa > pox f Px = a 

6 in hppt://www.nist.gov/wi 
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mqnd=mean(q(:,17)); 

Sqnda=stdiq(:;1/)) + 
qnds=(q(:,17)-mqnd) /sqnd; 

g3nd=mean (qnds.%3); 

g4nd=mean (qnds.%4); 

disp (‘Mean Std = direceroniess *) 
[mqnd sqnd] 

disp ('Skewness Kutosis - directionless') 
[g3nd g4nd] 
% ase 


(f'Directioniess wind speed:',' Mean=',intdstri{mand),' 
LE eS oe) |) 


a eae speed (moh}t} 


al_gumbel=pi/sqrt (6) /sqnd; 
u_gumbel=mqnd-0.577216/al_gumbel; 
eqamrn(g (2, 1) ) -lemax(gq (2,17) ) > 
fqq=al_gumbel* (qq-u_gumbel) ; 
Eqqg-aiecqumbel*exp (-fqq-exp({—-fqq).); 
figure 

hrseesi( ql); 30) 

title(['Directionless wind speed:',' Mean=',int2str(mqnd),' Std=',int2str(sqnd) ]) 
xlabel('Wincd speed {mph)}') 

ylabel (‘Histogram & Gumbel model') 

hold 

plot (qq, fqq) 

disp ("Gumbei al gumbel u_gumbei ') 

Le aloes u emer 


% Option i: Set sii im as, 373 


Shit cm mas > Li); 
shiftl=shift; 
qqmin=q(:,17)-shift; 
mqgqnd=mean (qqmin) ; 


iit 


sqqnd=std(qqmin) ; 

lamgl=mqqnd/sqqnd*2; 

kql=mqqnd*lamq1; 

qg=min:(q(*:,<b7).) “lemaxtq (ttt) 

qqm=qq-min(q(:,17)); 

fql=lamql*kql1* (qqm.* (kql-1) ) .*exp (-lamq1*qqm) /gamma (kq1) ; 
figure 

Histeest (a ea/g30) 


title(['Directionless wind speed:','. Mean=',int2str(mgnd),' Std=',int2str(sqnd) ]) 


xlabel ('Wind speed (moh) ') 


ylabel('Histogram & Shifted Gamma (Options 1/2: solid/dotted lines) 


hold 


Spiot (aq, faq) 


% Oisteseiny Sie aera con! 
‘ 3 momen 


skewness 


from kurtosis 


ae 4/g3nd*2 ; % 
% kqZ=2/ (g4nd/3-1); 
lamq2=sqrt (kq2) /sqnd; 
shift2=mqnd-kq2/lamq2; 


fom ‘aap eee e maa teas bpaktoan ah 3 
5 Set SERS Sloat Che CeO BES ER LOE Wee mich 


Ene. CONG LOnS Shi tiZ<min tats i ae 


Pete S eet Sineletater 
ShiftZ=shite; 
kq2=kql1; 
lamq2=lamq1; 
elseif shift2<0, 
shift2=0, 
lamq2=mgqnd/sqnd%*2; 
kq2=mqnd* lamq2; 
end, 
disp('Gamma (01) shee Ea Lamode') 
fShiit kql lamql1] 
disp('Gamma (02) shift k jambda') 
tshiii2 kq2 lamq2 ] 
qq2=qq-shift2; 
ae lamq2“kq2* (qq2.* (kq2-1)) .*exp (-lamq2*qq2) /gamma (kq2) ; 
ide ps 8 gh oe Tas 
=ctionless wind speed:',' Mean=',int2str {mond}, 


ce; z 
C 

Ses ee ae and) })} 

% 


2i{'Wind speed (mph) ‘3 
eai('Histogram & Shifted Gamma model') 


2HOLG 


%S WEIBULL DSITRIBUTI 


% (Method of moments) 
a a nl ce aa a eee 
% Estimated moments 


een ’ 
mw=mean (ww) ; 
Sw=std (ww); 
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m) 


wws= (ww-mw) /Sw; 
gw3=mean (wws.%3); 


% Mw Si 

a wee ee wee wee wee nee nes ~ - ~ 

% Skewness coe Snes 
o Po Sas Ste 3 oe arms 

% ch 19%, - Cnn, cma 

ee stare (Sk Hans (otis nses) Sn is enfin se Sin eS oe os 

OME Se, ee oe els Perna ey kB seth Ati wip lee Papa on 


dc=(cmax-cmin) /nc; 
ec=cminsdc: cmax; 
lc=length (cc) ; 

gl=gamma (1./cc+1) 
g2=gamma (2./cc+1) 
g3=gamma (3./cct+1) ; 
ae bee ROO le SO e ie ao oie 


Nie eS 


Cos 
a 


Geo 


ae 


cw= interpl (skew, CC; elles ‘spiline') 
noni gamma (1. Aero yee 

ggw2=gamma (2./cwt1) 
ggw3=gamma (3./cwt1) 
alw=sw/sqrt (ggw2-ggwl1%2) ; 
xiw=mw-alw*ggwl1; 

disp ('aipha G Sas L-RECORD! ') 
f[alw cw xiw] 


Neo 


No 


% Meso. ot and histogram for [-RECORD) 
VEXEwW 2 sD OF 

yw= (y-xiw) /alw; 
fw=(cw/alw) * (yw.* (cw-1)) .*exp (-yw.*cw) ; 

cdfw=exp (-yw.“cw) ; 

% Ligure 


SUPES oS” feeraye A ROUT a, 
ESE LWW, by SU} 


Ink eneedse?_? i oat r imam) 2 
GitA SMSO ° pS es Bo ° be (Magn } : 
uy 
} 
Vests & weeisi Hone 
Vistroaram for 
G See Oa, Oe 
5 eee CEaiey 
Pa a ea ee a ao aed Eat Se as aca han ar ia. ines A Hah VA can aR a a aA a VSS Co) Sb a a oS OR ahs esl 
& 


% REVERSE WEIBULL pdf and 


% histog Ler {RECORD) 


Nant! 


figure 
hist_est (-ww, 1, 30) 
hold 

piot (-y, fw) 

% titie('Reverse Weibuli fit and histrogram for RECORD') 
axis([min(-ww) max(-ww) O max(fw)+.01]) 

pees wind speed:',' Mean=',int2str(mqnd), Sta= pintestr(sqnd) |) 
xlabel ("Wind speed (m 

ylabel (‘Histogram & R 


feiouLl model!) 
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Sg tid RASRE 


Ze DART Tie, Colt 0s We ar ia OL KE fe yen geen es on Poy een Ym eee} 
% PARETO Di RIBUTION ea Dat GeO, na ey a ey 


% Construct, alwind sequence wish walues>=a\paréerco 


eg lab pieas, 
ti g (iL, 17)--at pareto, 
npar=npartl; 
qpar (npar) =q(i,17); 
end, 


5 Z Py: s\n) 

% Method of ts, based on 

om th plese es mrs Pe X? st = = — 3 a _ aaa Shae lGH 
% Buireme V eA DY H.. Case ON Hovomers eee, Ole. vv) 
2 ey a ean TR Nc a gy Tg ee ry a ar Co nS eRe SP Pp Gy PS eae Secale EN ey SG oa PE yes 
roy 


Gpar=qpar-a pareto, 

mpar=mean (qpar) ; 

spar=std(qpar); 

ky paretoi=((mper/spar)* 2-1) 727 

al_paretol=mpar®* ((mpar/spar) *2+1)/2; 

disp Sanpre-si72 k paretol aif parero) 
ee kKiparetoi altparetol] 


o Evie BSR el gr ns ee aes Pee Rene Hes Pe aero See 
% Probability weighted moments method, based on 

S 3t ewe Tae IE 3t Mx ey ‘or ee _- —~ 4 co ay die 

<S Extreme Vaiue ..., DY i MUS Sten Eo aeons Oe eae 

sic eat cet aa eat See ON Oe mc md Seta nse a ae psa a gs cit Sh San Gas She cs Sos a acs > an Eas OS Canes OAR a Ca ade ag A a CaN SO | coe oa CGE AC 
& 


qpars=sort (qpar); 

ppars=1:1l:npar; 

ppars=(ppars-0.35) /npar; 
tpars=sum((1l-ppars) .*qpars) /npar; 

k pareto2=(4*tpars-mpar) /(mpar-2*tpars) ; 
al_pareto2=2*mpar*tpars/ (mpar-2*tpars) ; 
disp('Sampie size e Vearevos ie 
[npar Kepareuol al _pareto2] 


& 
% astimate 

; a we a Pens ae cee ee 

“6 {from NIST site hppt://www.nist.gov/wind) 

¥ : 

I nr ee SY eye ec a Ss 


nparl=npar-1; 

mparl=sum(log(qpars (2:npar)+a_pareto) -log(qpars (1)+a_pareto) ) /nparl; 
mpar2=sum( (log (qpars (2:npar)+a_pareto) -—log(qpars((1)+a_pareto))=%2)/nparl; 
k pareto3=—(mpari+1-0° 5/ (l-mparl*2/mpanrz)); 


if eigcancs oe 

alt paretos=—appareto-~mpar Ll, 
else; 

al pareto3=4a-pareto*mpert aU +k pare.os) 
end, 
disp(‘'Sample size k pareto3 ai) parenos) 
[npar k pareto3 aly pareto3) 


*e _ =~] ~_ eee owe ee 
% DUERIN SS) ANE 

% moments 

% zy weighted moments 
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Pareto 3: De Haan method (NIST) 
figure 
HEStPES EA eal SU) 
title ([*Directioniess wind speed:*,'  Mean=",int2str(mqnd),' Std=',int2str(sqnd)]) 
xlabel (‘Wind speed (mph) ') 
ylabel ('Histogram & Pareto~i') 
LEK  paretoleo, 
xparl=a_pareto:l:al paretol/k paretolta_ pareto; 
xxpar=xparl-a _ pareto; 
£ paretol=(1-k_paretol*xxpar/al paretol) .*(1/k_paretol-1)/al_paretol; 
hold 
plot(xpari, £ paretol) 
elsei7) keparetol<0; 
xparl=a_pareto:1:max(qpar) ; 
xxpar=xparl-a_pareto; 
Peparetol=(1-k paretol*xxpar/al paretol) .~(1/k paretol—1) /al paretol; 
hold 
plot (xparl,f paretol) 
else, 
xparl=a pareto:;limax(qpar); 
xxpar=xparl-a_pareto; 
£f paretol=exp(-xxpar/al paretol)/al paretol; 


hold 
ploE(xXparl;£ paretol) 
end, 
figure 
hestrest(q 7730) 
title([’Directicniess wind speed:',' Mean=',int2str(mqnd),' Std=',int2str(sqnd) ]) 


xlabel ('Wind speed (mph)') 
Vlabel(*Histegram & Pareto~2') 
af Ke paretozed, 
xpar2=a_ pareto:l:al_pareto2/k pareto2+a pareto; 
xxXpar=xpar2-a pareto; 
Pepareve2-(1—k paretol *xxpar/alvpareto2) > (1/k paretoZ-1) /al paretoz; 
hold 
ploU(xpar2, = pareto2) 
elseitek=paretoz<0, 
RDate=appareto:1s>max(qpar); 
Rxpar=xparZ2-a pareto; 
fe Pace Lol— (it -kepareto2exxpar/ aly paretoZ yw (l/k pareto2—1)/aleparetoz; 
hold 
plor(xpar2, fi paretozZ) 
else, 
apase=alpareto:] :max(qpanr); 
xxXpar=xpar2-a_pareto; 
f pareto2=exp (-xxpar/al_pareto2)/al_pareto2; 
hold 
plot (xpar2,f pareto2) 


figure 

histeest(q;,.17;30) 

title(['Directioniess wind speed:',' Mean=',int2str(mqnd),' Std=',int2str(sqnd) ]) 
xlabel ('Wind speed (mph) ') 

ylabel (‘Histogram & Pareto~3') 


ibs 


Lf k paretoi>0, 
xpar3=a_pareto:1l:al_ pareto3/k pareto3+a_pareto; 
xxpar=xpar3-a pareto; 
f pareto3=(1-k_pareto3*xxpar/al pareto3) .*(1/k_pareto3-1) /al_pareto3; 
hold 
plot (xpar3,f pareto3) 
elseif ko pareto3<0;, 
xpar3=a_pareto:1:max(qpar) ; 
xxXpar=xpar3-a_pareto; 
f pareto3=(1-k_pareto3*xxpar/al_pareto3) .*(1/k_pareto3-1) /al_pareto3; 
hold 
plot (xpan37f paretes) 
else, 
xpar3=a_pareto:1:max(qpar) ; 
xxpar=xpar3-a_pareto; 
f pareto3=exp (-xxpar/al_ pareto3)/al_pareto3; 
hold 
plot (xpar3,f pareto3) 


MONTE CARLO GENER? ee OF RURRICANE HAZARD 
% AY THE PM 


MT MITES a ee Cegees res 3 exe + ieee 
Al, TIMES (No Gistinction is meade between 


% the hurrican and non-hurricane seasons) 

rand('seed*,nseed) 

time=0; 

ktime=0; 

while .time<=nyr, 
ktime=ktimet+l; 
time=time-log(rand(1,1))/nu; 
thr (ktime) =time; 

end, 

nhurr=ktime-1; 

ilabiare—elare(( Ike ilalivieie)) 


Reverse 


% ShiftedsGamna,  Opbicnemimand 2 


x enue mc=shiftl+icdf ('Gamma',uu,kql,1/lamql1) ; 
x shg2_mc=shift2+icdf ('Gamma',uu,kq2,1/lamqZ2) ; 
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Pe keparecol==(, 
= pearl eme=arpareto-al sparetol* tog (uu) 
else, 
Beparieme—aupareét 04a IWparetol*(i- (uu)Pekeparetol)/ k paretob; 


se 


end, 
tiekepareto2==0, 

Mapabe mc=4spareto—al |/pareto2* log (uu) 
else, 

sapere wC=aeparTetcoral pareto2”~(1=(au)iaek paretoZ)/ k-parete2; 


‘Ne 


end, 


Peekeparcioo==0, 
Mepans mMc=aypareto ral spareto3* log (uu) 
else, 


Sepaco mC=anparecotal, pareto3s* (1— (um)r *keparetos)y/ k paretos; 


“Me 


AVERAGE 


eS 
¥ iC Gig 


& PE eRTODS 3S ak 


THE ; 
% ON EMPTREICAL DISTRIBUT 
Cate 50) 2 iL & IMOOO Os 

nutau=nu*tau; 

56 nd=q (:,17)% 


x stairs ixg 3 (xg = eee 


Generated Gumbel, shifted Gemma, Reverse Weibull, 
é and Pareto records 
p 
z Gain 
% Guripe 
Cc 


ledf-nd,xq nd, loc, upcl. = ecdf (x gumbel mc)’; 
x gumbel=interp] (cdf_nd,xq_nd,1-1./nutau, ‘spline"); 


eae no, “pind, Loc, upc} = ecdi (x shglemc) ; 
x shgaml=interpl (cdf_nd,xq_nd,1-1./nutau, 'spiine') 


“ee 


LF 


[edi nd; xqund, loc, upc) ="ecai{xtsngzimc);, 

x_shgam2= eae ee nd, xqend; L=1 .7nutau, ter sere) 
% xe Senge i=shifti+s ea emma IE Genes: Lge i 
: "Gamma 


means a 
[cdi nd, xq nd, boc,upe |p = ecdé (xarwemc)i, 
x rweinterpl (ed£ nd; xq nd, t—-1- f/outau, Sspiane" )e 


% iaiw Cw xe Wi 

& 3 Wes WD py ay Pn ena ee BET acer Meee ee 

te HEOee { "wo 7 oe eae Ew le 

(ou 

res gt wees a on = ~ es t t 

% shgam, tau, x _shgar 2, ew Ay 

% vyerage return period ect eb 

% REVS SSCSee 

rom aa $ BT ac a EE ote PON 

% ~solid Weibuli-~dotted line’) 
in ate eaee es 2 

é Parezeo.  Opirens 8,2, ane 3 

a 


[ede nd,xq_nd,loc, aiten ="ecdt (x epanl me)y 
x parl=interpl (cdf _nd,xq_ nd,1-1./nutau, ‘spl: 
[cdf nd,xq nd, loc,upe] = ecdfi(x*par2z*mc); 

x par2=interpl (cdf nd;7xq nd, l-17/nutau;. sone %) 
[cdfind, xq nd, loc, upc |) =seed h(x parsameyT; 
x Poets meatier = ndyescq end lS layinuta, az 


Lig) 
b 

} 3 
2) 
® 
~~ 
“eo 


“ee 


spline’); 
(2. / DUE Al) 


Lomita). % 


2 Frat oe ear So \ © 
boa / (aE 


ge oe ae de 


A LD a eed pain ye een te en OK Taher ge arene Bond eae eet ANY 
S Wei aS Pari e area Goes iOee {i EB an fh lok cbete eb te 
tg 
eee tle Eee wat i Pap ens ht sew GRRE) SR eR Ff i3 { Rae aie 
Une Oke ctiy Ors SOs wee ss 2 Avr Sabiy: © 
Be es (Lae Gone Bice eT we Fig / Poa S| ENGR RAR aA SUN Sf a, Sane 
SUOrTrat Paetiocegs Sek. (MWe yy Neat SiO fi cece ees ey, 
e = 
ve = = seo stn se se wesees see sses ps so 


Yeers 
4 


WIND Be a 


EVERED BY SE 
RSE WEBULL (x ant: and 


Hh? of ale 


% = ESTIMATE BASED ON 
a. 
a cacao Nl tl lt a ct oa 


figure 

ploti(tau xrgumbel et .s | causes wy yale pars, call, CMp yeas) 
xlabel('Average return period (years)') 

ylabel (‘Wind speed (mph)'") 

title ('Gumbel (dotted); Reverse Weibull (dashed); Pareto-DEHAAN (solid); 
Empirical (dashdot) ") 


figure 


plot (tau,x_shgaml, , tau, x: parl;, ==", tan, xi par2;, P24; taupxk par3;itau, xvenp, ee 
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ae 


xlabel('Average return period (years) ') 

ylabel (‘Wind speed (mph)') 

title('Gammali(o); Reverse Weibuli(dashdot); Paretoi;2;3 
Empirical (dashdot) '") 


a 


a 


. elt BS, 
‘dashed maAttred«=aniad) = 
aS MG CHO eee Ges SIO a G3 


& 2D . + - : pig ae 7 

G Reverse Weibul WEED se SWS 
oo x peut > 2 ete a ah a ede oA EB en SSID dca ee ee 
oe 


xX Irws_mc=—-xiw-icdf ('wbli',uu,alw,cws) ; 
Pcdtend, xq nd, loc,;upe] = ecdf (x xrws me)’; 

x rws=interpl (cdf nd,xq _nd,1-1./nutau, *spline') 
figure 

BP loe(tau;, x rw, -—', tau, x rws,,tau,x emp,’ -..*) 
xlabel (‘Average return period (years)'") 
ylabel('Wind speed (mph)') 


“Ne 


9 


title('Reverse Weibull: estimated(dashed) & imposed(solid) tail parameter; 
Empirical {dasndot)*) 

o. 

2. 


Pthurr,x gumbel mc, x_ 
% nd mc (.5,30, LQ00, 
reto = 30; 40 


year, > S parl hiG, \ peas Mey xe ers mci=hurer 
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